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^ (— I We present a numerical study of the time-dependent SN equations 

in three dimensions with three kinds of symmetry: spherically sym- 
metric, axially symmetric and translationally symmetric. We find that 
the solutions show a balance between the dispersive tendencies of the 
Schodinger equation and the gravitional attraction from the Poisson 
^ ■ equation. Only the ground state is stable, and lumps of probability 

^ . attract each other gravitationally before dispersing. 

1 Introduction 

Ths is the second in a series of papers presenting a numerical study of the 
Schrodinger-Newton (or SN) equations. In the first ([|J) we reviewed known 
results on the SN equations and then analysed the linear stability of the 
spherically-symmetric stationary states. We found that the ground state 
is linearly stable, in that all the eigenvalues for linear perturbations were 
purely imaginary, while all the higher states are unstable. The (n + l)-th 
state, or equivalently the n-th excited state, has n quadruples (A, —A, A, —A) 
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of complex eigenvalues with nonzero real part. We now want to go further 
in two directions. On the one hand we shall consider the full nonlinear 
evolution of the SN equations and on the other we shall move away from 
spherical symmetry, allowing for two space variables. 

In the evolution of spherically-symmetric data we shall find that the pic- 
ture from linear theory is confirmed: the ground state is stable but slight 
perturbations of the higher states decay, with probability density escaping 
off the grid to leave a 'nugget' consisting of the ground state rescaled to have 
total probability less than one and situated at the origin. For more general 
spherically-symmetric data the picture is confirmed: solutions typically dis- 
perse leaving a nugget of rescaled ground state. (By a rescaled ground state 
we mean the following: if ipo{r) is the spatial part of the wave function for 
the ground state, normalised to have total probability equal to one and with 
conserved energy So in the terminology of ||, then p 2 ipo(pr) with p a posi- 
tive constant, is a stationary solution with conserved energy p 3 £ and total 
probability p.) 

We can introduce a new space variable in two ways. We may consider 
axially-symmetric solutions in 3-dimensions. In this case there are new 3- 
dimensional stationary states, somewhat analogous to the axially-symmetric 
solutions of the hydrogen atom, and with nonzero expectation for the angular 
momentum. These turn out to be unstable under evolution, as one would 
expect. For the solution most analogous to a pure dipole, the evolution 
clearly shows the two 'lumps' of probability falling from rest into each other 
and evolving towards the ground state, with some scatter. 

Alternatively we may consider translation-invariant solutions in 3-dimensions 
or equivalently the SN equations in 2 + 1- dimensions. (Of course the solutions 
are not normalisable as 3 + 1-dimensional solutions but this is still an inter- 
esting problem.) In this case we can find rigidly- rotating stationary states 
like two lumps of probability rotating around each other. When these are 
evolved with the time- dependent SN equations they prove to be unstable, 
so that it is possible to arrange for two lumps of probability to be in orbit 
around each other at least for a while before they merge. 

In summary, the picture that we find is of a system with dispersive ten- 
dencies because of the Schrodinger equation and attractive or concentrating 
tendencies from the gravitational attraction. There appear to be infinitely 
many stationary states, all of them unstable except for the ground state. 
General data evolve to leave some residual probability in a rescaled ground 
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state with the rest of the probability dispersing. Lumps of probability in the 
initial data can attract each other and even orbit each other but eventually 
the dispersive tendencies win. 



The plan of the paper is as follows. We shall end this Introduction with 
an analytic result bounding the residual probability left on the grid at late 
times. In Section 2 we introduce the numerical method which we shall use to 
evolve the spherically-symmetric SN equations. The results of this evolution 
are presented in Section 3. In Section 4 we present the results of numerically 
evolving 3-dimensional axisymmetric solutions and in Section 5 we present 
the results of numerically evolving the 2-dimensional SN equations. 

For the analytic calculation then, suppose that the picture presented 
above holds in general so that arbitrary initial data evolve to give a scat- 
tering solution which disperses to infinity and leaves a rescaled ground state. 
Then we can obtain a bound on the probability remaining in the ground 
state provided the initial energy is negative. For suppose the initial value of 
the conserved energy is Si then since this is conserved, we can calculate at 
late times to find 



where 8$ is the energy in the scattering solution, which we suppose to be 
positive, So is the (negative) energy of the ground state and p is the proba- 
bility left in the ground state. (The ground state is strongly peaked at the 
origin so that the cross-term in Sj is zero.) Taking account of all the signs, 
if Sj < we can rearrange this to read 



2 Numerical methods for the spherically-symmetric 
evolution 



S I = Ss+p 3 S 




(1) 



which is the desired bound. 



Recall from j| that the SN equations are the system 



(2) 
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In this section and the next, we are concerned only with the spherically- 
symmetric case and in this case (|2]) can be simplified as 



.du d 2 u 
dt dr 2 
d 2 (rd>) \u\ 2 



dr 2 r 



+ <fm, (3) 

(4) 



where u = rip. 

To solve this system, we shall use a Crank- Nicholson method for the time- 
evolution of the Schrodinger equation and the spectral method of for the 
Poisson equation. Schematically, this takes the form 

n+l _ n 

2i - = -D 2 (u n+1 ) - D 2 (u n ) 

ot 

+0"+ V +1 + <f) n u n , (5) 

where, with the spectral method of 0, D 2 is the second- derivative matrix 
for Chebyshev polynomials. The Crank-Nicholson method is second-order 
accurate in the time and preserves the normalisation of the wave function. 

Now we need an iterative method to find n+1 . We write 0£ +1 and u^ +1 
for the iterates and make the initial choice 0o +1 = 0™, then use the spectral 
method to solve fl5|) for Uq +1 , use this in (|J) to improve 0q +1 to 0i +1 , and 
solve (d) to find This cycle can be repeated until the desired degree of 

convergence is reached. 

The boundary conditions are that u and rcf> should vanish at the origin 
and the outer edge r = L of the grid. We expect there to be an outgoing 
flux of probability and we need to prevent it from reflecting back off of the 
boundary. For this we use the 'sponges' of Bernstein et al [Q. The idea 
is to replace the i on the left-hand-side in (^) by i + s(r) where s(r) is a 
positive function large (of order one) near the boundary but small close in. 
We typically use s(r) = ae b ^ r ~ L ^ for constants a and b. By this means the 
Schrodinger equation is converted to a heat equation near the boundary and 
the outgoing flux is absorbed. 

We can perform various checks on the method outlined here. To check 
that the sponges are working as desired, we temporarily set = and con- 
sider the numerical evolution corresponding to the following explicit spherically- 
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symmetric solution of the zero-potential Schrodinger equation: 



rip = u 



[exp( 



r — vt — a) 2 ivr 



ivH 



(a 2 + 2irf 



2(a 2 + 2it) + 2 



4 

iv 2 t 



exp( 



(r + vt + a) 2 zt>r 



)] 



(6) 



2(a 2 + 2it) 2 
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where C is a normalisation constant (which can be found explicitly). This 
solution is a spherically-symmetric bump with a Gaussian profile centred 
initially at r = a and moving radially with velocity v. It both disperses and 
moves off the grid. Computing it is a test for the Schrodinger evolution, and 
having it move smoothly off the grid is a test for the sponges. The initial 
data for the solution (Jg) will be evolved in the full evolution in Section 3. 

Next, we can include a fixed in (|3|), dropping (f|) and jointly test the 
Schrodinger evolution and the sponges. Finally, we can test convergence of 
the method by varying the time step and/or the number of Chebyshev points 
chosen in the spectral method. 

3 Results in the spherically-symmetric case 

Having tested the method and the sponges satisfactorily, we evolve the 
ground state, that is we use as initial data the stationary-state of lowest 
energy which is known from |J or Q. Since this is a stationary state it 
should evolve with a factor e~ lEt and this is just what we find (see figure |1|: 
the phase is linear in time). However if this evolution is continued for a 
long enough time, numerical errors accumulate and the solution drifts away 
from the ground state. By Fourier analysing the solution one finds lines in 
the power spectrum at (approximately) the normal frequencies for a pertur- 
bation about the ground state already found by linear perturbation theory 
in Q. Thus we really are getting the stable but perturbed ground state. 
We have also evolved from data which are the ground state with a small 
perturbation introduced by hand and the picture is the same. 

Next we evolve the second spherically-symmetric stationary state taken 
from Q. We should expect this also to evolve as a stationary state at least 
for a while until numerical errors accumulate. What we obtain is shown 
in figure |[ This shows \rip\, so that the second state, which has a zero, 
appears as bimodal. The evolution proceeds as a stationary-state to about 
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time 



Figure 1: The phase angle of the ground state as a function of time. 

t = 2000 and then there is a sudden change to a unimodal solution with some 
oscillation. This we claim is a perturbed and rescaled ground state solution, 
with total probability less than one. 

To support this claim, we have again calculated the power spectrum and 
find lines at (approximately) the normal frequencies of perturbations about 
the ground state (to find agreement we need to rescale these frequencies with 
the factor corresponding to the rescaling of the residual ground state). We 
can also use the check suggested by equation ([!]). That is, we calculate both 
the probability p and the action or conserved energy £ remaining on the grid 
at time t and compute the bound on p provided by (|l|). If the argument 
leading to this bound is correct then these two probabilities should converge 
on each other and as we see in figure they do. Again, we can evolve from the 
second state with a small perurbation introduced by hand rather than waiting 
for the numerical errors to accumulate when the collapse to the ground state 
is more rapid, or alternatively we can evolve from the second state but with 
a more stringent tolerance on the iteration determining <f>, when the collapse 
takes longer. 

Moving on, the decay of the third state is shown in figure The initial 
state this time has three maxima and evolves for a while as a stationary 
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time 200 r 

Figure 2: Evolution of second state. 
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Figure 3: Decay of the second state: the two measures of probability con- 
verging 
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Figure 4: Evolution of third state. 

state before collapsing to a rescaled ground state with the emission of some 
probability to infinity. We have plotted the figure corresponding to figure § 
and we again find convergence. 

The picture that emerges from these calculations is of nonlinear insta- 
bility of all states after the ground state. Each higher state decays to the 
ground state either by the accumulation of numerical errors or because of an 
explicitly included perturbation. The end result is a noisy rescaled ground 
state, noisy because the eigenvalues for perturbation about the ground state 
are all imaginary. 

The other spherically symmetric evolutions which we have calculated are 
with the initial data furnished by equation @ with t = 0. This is a Gaussian 
bump centered at r = a with width a and moving with velocity v. The 
conserved energy of this data rises with v (because the kinetic energy rises) 
and with a (because the particle is more localised). In all cases the solution 
disperses leaving a rescaled ground state which we can characterise by the 
probability residing in it. In figures [5], ^| and [7] we plot the residual probability 
against v, a and a at different times. 
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Figure 7: Evolution of Gaussians with varying a: v = 0, a = 50. 



From figure [5] we see that the residual probability is greatest if the lump 
is released from rest, but it is easier to trap probability if the initial velocity 
is ingoing than if it is outgoing. From figure || we see that for lumps released 
from rest, the residual probability is greater if the lump starts closer in, or 
in other words is more gravitationally bound. From figure [7| we see that 
reducing a, which raises the energy, leads to more dispersion. 

A key motivation for this set of calculations was the desire to check con- 
vergence of the method, and we have done this in a variety of ways. First for 
the time-step, with a Crank-Nicholson method we expect to have quadratic 
convergence, and this can be investigated with a Richardson quotient. We 
suppose that for some variable of interest the calculated value Oh and actual 
values I are related by 

O h ~I + Ah\ 

where the time step size is h, and k is the order of the error, then we can 
calculate the Richardson quotient: 

— Oh 3 
Oh 2 - h3 ' 

where hi,h,2 and h% are three different values of the step size. With k known, 
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this will be a simple function of the hi and so will provide a check on k. Next 
for convergence in space we can repeat the calculation with different numbers 
N of Chebyshev points. For both checks, the results converge as required. 

With this calculation we confirm the picture of the SN evolution which 
we have been claiming: the Schrodinger equation tends to disperse probabil- 
ity but the Newtonian gravitational attraction holds it together; all states 
after the ground state are unstable and the evolution in general leads to a 
dispersion of probability to infinity, leaving some residual probability in a 
rescaled ground state. 



4 The axially- symmetric SN equations 

In this chapter, we solve the SN equations for an axially- symmetric system 
in 3 dimensions. The wave-function is now a function of polar coordinates r 
and 9 but is independent of polar angle. We shall first find stationary solu- 
tions. These typically have nonzero total angular momentum and include a 
dipole-like solution which appears to minimise the energy among wave func- 
tions which are odd in (the usual) z. We then consider the time-dependent 
problem. In particular we evolve the dipole-like state, and it turns out to 
be nonlinearly unstable - the two regions of probability density attract each 
other and fall together leaving a multiple of the ground state as the evolutions 
in Section 2 did. 

With u = rip the system of equations @ becomes: 

-u„ 1 (sin 9ug) e + 4>u, (7) 

r z sm 9 

(r<p) rr H ^— (sin 9(f) 9 ) e . (8) 

r sm 9 

For stationary solutions, the left- hand- side of (^) is replaced by Eu. Bound- 
ary conditions are that u = at r = and r = L, = at r = L and finite 
at r = 0, and ug = = <pg at 9 = and 9 = it. 

To find stationary solutions we proceed as follows: 

1. take as an initial guess for the potential <fi = -. -; 

(1 + r) 

2. using this potential solve the time- independent Schrodinger equation; 



iu t = 
u\ 2 
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3. select an eigenfunction in such a way that the procedure will converge 
to a stationary state (this needs trial and error); 

4. calculate the potential due to the chosen eigenfunction; 

5. take the new potential to be the one obtained from step [| above but 
symmetrised, since we require that should be symmetric around the 

7T 

6 = — (otherwise numerical errors can cause the wave-function to move 



along the axis); 

6. now provided that the new potential does not differ from the previous 
potential by some fixed tolerance in the norm of the difference, stop, 
otherwise continue from step 2. 

In step 2 we solve the eigenvalue problem by a 2-dimensional spectral 
method, using Chebyshev differentiation in the directions of r and 9. In step 
3, the eigenfunctions at the first iteration are labelled by the usual l,m and 
n quantum numbers, though with m = for axisymmetry. The idea is to 
choose one and run through the iteration, hoping for convergence. When the 
method converges, we do obtain a stationary state but we do not arrive at a 
one-to-one correspondence between solutions of the starting linear problem 
and the final nonlinear problem. 

For the energy eigenvalue E (which is not the conserved energy) we have 
the formula 



while the total angular momentum J 2 , after integration by parts, is : 



We present in Table [1| the first few stationary states of the axially sym- 
metric solution of the SN equations ordered by their energy and named for 
convenience axil to axi8. 

In this table, axil, axi3 and axi8 are spherically-symmetric solutions turn- 
ing up again (as they should). Axi2, shown as a surface in figure ^]and as a 
contour plot in figure is very much like a dipole solution and appears to 
be the solution minimising the energy among wavefunctions odd in z. If so 





(9) 




(10) 



Both integrals are over R 3 . 
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Figure 8: The dipole-like state, axi2. 
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Figure 10: Contour plot of the state axi4. 




Figure 11: Contour plot of the state axi5. 
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Figure 13: Contour plot of the state axi7. 
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Energy 


J 


name 


-0.1592 


zero 


axil 


-0.0599 


5.1853 


axi2 


-0.0358 


0.002 


axi3 


-0.0292 


2.3548 


axi4 


-0.0263 


17.155 


axi5 


-0.0208 


3.1178 


axi6 


-0.0162 


5.2053 


axi7 


-0.0115 


1.9E-6 


axi8 



Table 1: The first few axially-symmetric stationary states, 
it was found already in || . Contour plots of axi4-axi7 are given in figures [H] 
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To solve the time-dependent axisymmetric SN-equations, we shall use an 
alternating direction implicit (or ADI) method (see e.g.|J, @ or 0) We split 
the Laplacian in ([7|) to write it as 

u = i{L x + L 2 - 4>)u (11) 

where 

1 d 2 

Ll = rd? (12) 

U = ;[^ + cot^]. (13) 

Next we introduce new variables S and T and write this in the formally 
equivalent form: 

ex P(--y£i)S(» = exp(yL 2 )w(t) 

exp(-^L 2 )T(t) = exp(yL 1 )S(t) (14) 
%h %h 

ex P(y + h ) = ex P(-y^) T (*) 

where we suppress the dependence on r and 0. To obtain a discrete form of 
(|T4l ) we linearise in the time-step h to find 
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(1 + 

(1 - |^)T'< 



'it 



m 



(15) 




which we solve at each instant by a Peaceman-Rachford ADI iteration |T|. 
This means we discretise it as 



and iterate with 0° = 0. Here k labels the iteration (at a fixed time) and p 
is a small constant chosen so that the iteration converges. 

The boundary conditions are that u and vanish at the outer bound- 
ary and so we need sponges as in Section 2 to prevent waves of probability 
reflecting back. 

As an example, we take as initial data the dipole state of figure |[ The 
evolution of | -0 1 for this is shown in figure [TJ. (Note that -0 is initially real 
but must become complex before again becoming approximately real in the 
remote future. The initial data for ip is odd as a function of z but it ends up 
approximately even.) 

Following the lessons learned in Section 3, we should expect the stationary 
state to be unstable. We can compute the evolution for short times to see that 
initially it remains a stationary state, that is that the only change is a phase 
growing linearly with time. However when numerical errors have built up the 
state becomes unstable and the two concentrations of probability density fall 
into each other. Probability leaves the grid, as does angular momentum, and 
we are left with a multiple of the ground state. We can check convergence of 
the method by calculating a Richardson quotient with three different time- 
steps (to find that convergence is now linear in time), and by varying the 
number of Chebyshev points. 

What we learn from this calculation is that, as well as dispersion which 
is what we mostly saw in Section 3, the solutions of the SN equation show 



(Li + p)4> k+ i 



(L 2 - p)4>k + ~^\u 
(Li - p)(j) k + -^\u 



(16) 



(L 2 + p)(pk+i 
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t = 1400 







Figure 14: Evolution of the dipole. 



gravitational attraction: lumps of probability density released from rest fall 
into each other. In the next section we shall see that, at least for the two- 
dimensional SN equations, lumps of probability can orbit each other. 



5 The two-Dimensional SN equations 

In this section, we shall consider the SN equations in a plane, that is in 
Cartesian coordinates x,y. We shall find a dipole-like stationary solution, 
and some solutions which are like rigidly rotating dipoles. These rigidly 
rotating solutions are unstable however and will merge, radiating angular 
momentum. 

The SN equations in this case are 



I1p t = -1p XX - 4>yy + (j)i> (17) 

'x.r i 4*yy 



+ 4>yy = M 2 



We use the ADI scheme as in equation (|T^) but with the understanding that 
now 

In = — 

dx 2 
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dy 2 



For the potential we use the counterpart of (jig ) with the same understanding. 
The boundary conditions are that ip and <fi vanish at the edges of the grid, 
which is now a large square. We still need sponges and we do not want them 
to have corners so we take functions like s(x,y) = min[l, e^V^ 2 ^ 2 ) -20 )]. 
We can test the efficacy of the sponges by evolving moving two-dimensional 
Gaussians to see that they propagate off the grid, as they do. 

Once we have the code, we can look for solutions corresponding to those 
found already in 3-dimensions. In particular we find a stable ground state 
and a dipole-like solution which is stable for a while before decaying to the 
ground state. We can seek solutions with no counterpart among those found 
already by making an ansatz of rigid rotation. That is we look for solutions 
which in polar coordinates r, 6 take the form: 

if,(r, 6, t) = e- iEt ip(r, 9 + ut) (18) 

for (real) constants E and cu. (Solutions of a similar nature in 3-dimensions 
would depend on all three spatial coordinates which is why we haven't seen 
them so far.) 

To separate the SN equations for a solution like QlSD , we go into the 
rotating frame with Cartesian coordinates X, Y given by 

X = x cos ujt + y sin ut, 
Y = —x sin ut + y cos ut, 

and the time-dependence separates off to leave the equations as 

- V 2 ^ + H - iu(Yi> x - XiJjy) = Eip (19) 

V 2 = |VI 2 , 

To solve ATP]), we begin with a small value of u and the dipole-like solution. 
Iteration leads to a solution, and we can study its change with increasing 
uj. The wave-function i/j is necessarily complex and we display the real and 



imaginary parts of such a solution, with u = 0.005, in figures [L5| and [16. 

The solution just found is a stationary state but we must expect it to be 
unstable. If we use it as initial data for the time-dependent problem then we 
find as before that it evolves as a stationary state for some time (see figure [l 
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Figure 15: Real part of spinning solution, uj = 0.005. 
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Figure 16: Imaginary part of spinning solution, uj = 0.005. 
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t = 



t = 200 




before numerical errors build up and it becomes unstable (figure |18|). The 
two lumps of probability orbit each other about four times before collapsing 
into a single lump in about one orbital period. Probability and angular 
momentum are radiated off of the grid. 

What we have found by this calculation is a confirmation of the picture 
of lumps of probability interacting gravitationally with each other. Here the 
lumps are in orbit around each other for a while before becoming unstable 
and collapsing into a single lump. By earlier work, we must then expect the 
probability to disperse leaving a rescaled ground state, as it does. 
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t = 5200 



t = 6400 




Figure 18: Later evolution of spinning solution: collapse. 
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